library(automap)
library(gstat)
library(dplyr)

#   automatically fit variogram with automap  
# set the formula of ordinary krigin
formula.volume_rate<-'volume_rate~1'%>%as.formula
valid.volume_rate<-ha_info.sp[!is.na(ha_info.sp$volume_rate),]%>%
    spTransform(CRSobj = crs(poi.diversity))%>%remove.duplicates
# zerodist(valid.volume_rate) make sure there is no duplicates
# random sample test set data
test.inx<-sample(1:nrow(valid.volume_rate),nrow(valid.volume_rate)*0.2)
training_set.volume_rate<-valid.volume_rate[-test.inx,]
test_set.volume_rate<-valid.volume_rate[test.inx,]
mm<-autofitVariogram(formula.volume_rate, training_set.volume_rate)
volume_rate.krige.fit <- gstat(NULL, formula =formula.volume_rate,
             data =training_set.volume_rate, model=mm$var_model)
# cv.volume_rate<-gstat.cv(volume_rate.krige.fit,nfold = 3)
# volume_rate.krige<-interpolate(poi.diversity, gUK, xyOnly=FALSE)
# compare.cv(list(krige.cv_output = cv.volume_rate))
# t1<-system.time(valid.volume_rate.pre<-predict(volume_rate.krige.fit,
                               # newdata=test_set.volume_rate[1:10,]))

t2<-system.time(valid.volume_rate.pre.2<-predict.gstat.par(volume_rate.krige.fit,
                               s=test_set.volume_rate))
# cat(identical(valid.volume_rate.pre,valid.volume_rate.pre.2))
# if(!require(hydroGOF)) install.packages("hydroGOF")
# library(hydroGOF)
# gof(sim=valid.volume_rate.pre.2$var1.pred,
#               obs=test_set.volume_rate$volume_rate)
# 
# library(intamap)
# training_set.volume_rate.small<-training_set.volume_rate[
#     sample(1:nrow(training_set.volume_rate),1000),]
# training_set.volume_rate.small$value<-
#     training_set.volume_rate.small$volume_rate
# anisPar <- estimateAnisotropy(training_set.volume_rate.small)
# print(anisPar)  
# rotatedObs <- rotateAnisotropicData(training_set.volume_rate.small,
#                                     anisPar)
# newAnisPar <- estimateAnisotropy(rotatedObs)
# print(newAnisPar)  
# intamap.obj <- createIntamapObject(
#     observations = training_set.volume_rate,
#     formulaString=as.formula(volume_rate~1),
#     predictionLocations = meuse.grid,
#     class = "idw"
# )
# estimateParameters.automap()